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Foreword 


The National Aeronautics and Space Administration and the Atomic Energy 
Commission have established a Technology Utilization Program for the dissemina- 
tion of information on technical developments which have potential utility outside 
the aerospace community. By encouraging multiple application of the results of 
their research and development, NASA and AEC earn for the public an increased 
return on the investment in aerospace research and development programs. 

This publication, part of a series to provide such technical information, is 
intended for mathematicians, engineers, physicists, and other persons with training, 
knowledge, and interest in mathematics. A small part of the large field of theoretical 
and applied mathematics is covered. Section 1 primarily concerns mathematical 
applications for the resolution of problems encountered in numerous industries. 
Section 2 concerns computer programs which can be applied to resolve computa- 
tional problems swiftly and accurately. 

The usefulness of this material is indicated by its applications to aircraft, 
automobiles, beams, cylinders, structural elements, oil fields, geophysics, missiles, 
plates, solid rocket propellants, mobile nuclear reactors, rods, solutes, solvents, 
space vehicles, spheres, pressure vessels, and ultrasonics. 

Additional technical information on individual devices and techniques can be 
requested by cifcling the appropriate number on the Reader’s Service Card included 
in this compilation. 

Unless otherwise stated, NASA and AEC contemplate no patent action on 
the technology described. 

We appreciate comments by readers and welcome hearing about the relevance 
and utility of the information in this compilation. 

Ronald J. Philips, Director 

Technology Utilization Office 

National Aeronautics and Space Administration 
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Section 1. Mathematical Applications 


CURVE FIT TECHNIQUE WITH FIRST DERIVATIVE CONTINUITY 
USING A MINIMUM NUMBER OF ANALYTICAL EXPRESSIONS 


This curve fit technique is a chain system of 
connected segments or arcs of the total curve, 
where this curve is single valued at all points with 
continuous first derivatives at these points. The 
analytic expression, y = A /( 1 -Bx) +Cx 2 + Dx +E, 
closely represents each segment. 

The coefficients A, B, C, D, and E are deter- 
mined by three points on each arc or segment and 
the two slopes of the first and third points, where 
these two points represent the ends of the segment 
or analytical expression. Determinations of the 
coefficients are given in a detailed development 
of the equations. 

The technique requires less analytical expres- 


sions with more accuracy than other currently 
used techniques; it also offers a simple, more 
direct method of solution of the coefficients from 
the curve fit for three points and the adjacent two 
slopes of the curve, and elimination of the compli- 
cations of wave forms, or excursions, of the curve 
between adjacent sets of points. 

Source: T. M. Dannback of 
The Boeing Co. 
under contract to 
Marshall Space Flight Center 
(MFS-14735) 

Circle 1 on Reader s Service Card . 


COMPUTATIONAL PROCEDURE FOR A MINIMUM-MAXIMUM 
POLYNOMIAL CURVE FIT ROUTINE 


Some form of least-squares curve fit is com- 
monly employed to fit a span of data or to ap- 
proximate a mathematical function for use in a 
computer. If the maximum error resulting from 
the approximation is of interest, however, this pro- 
cedure is not optimum. The maximum error may 
be reduced 10 to 20% by iteration on the curve 
fit. 

A mathematical procedure for improving fitting 
between a set of data points and a polynomial 
uses a minimum-maximum fit, i.e., one in which 
all the maximum errors are of the same size. This 
procedure uses a least-squares iteration; the de- 
pendent variable is the deviation of each local 
maximum (including end points) from the average 
value of the local maxima; the independent vari- 
able is the array of polynomial coefficients. 

An initial least-squares fit is assumed to be 
available. The error curve is computed. For a 
polynomial approximation of degree n, an error 
curve with n 4- 1 values of zero and n + 2 local 


maxima (including the endpoints) normally re- 
sults. The “average” maximum error (either the 
arithmetic average of magnitudes or root-mean- 
square) and the deviations from this average for 
each local maximum are determined. A correction 
to the polynomial coefficients involving minimi- 
zation of the extremes of the error curve is com- 
puted by a weighted least-squares method. The 
above steps from computation of the error curve 
are repeated until the deviations of the errors from 
the average are within some specified tolerance. 

The program is designed to remove extraneous 
local maxima and drive the curve fit to the desired 
solution. Isolated cases in which the endpoint 
is not a local maximum may also be treated. The 
logical flow to delete extraneous local maxima is 
believed to be novel. 

A method is provided for the best polynomial 
approximation in the sense that the error curve 
has the lowest maximum error; this is obtained 
at a loss of 1 to 2% greater root-mean-square 
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error. The maximum error is reduced 10 to 20%. 
There is a potential wide use for this method in 
computer science to evaluate complex geometric 
surfaces, such as the design of automobile bodies 
and aircraft fuselage and wing configuration. 


Source: J. B. Clifford, Jr. of 

TRW Systems Group 
under contract to 
Manned Spacecraft Center 
(MSC-13257) 

Circle 2 on Reader s Service Card. 


RATIONAL CHEBYSHEV APPROXIMATIONS FOR FERMI-DIRAC 
INTEGRALS OF ORDERS -1/2, 1/2, and 3/2 


The complete Fermi-Dirac integrals are usually 
defined by (x) = ./ 0 °° — — ,k>-l. These 

e ( l ~ x ) f~i 

integrals appear in a variety of applications sub- 
ject to Fermi-Dirac statistics, for example, in the 
theory of semiconductors. The most frequently 
used functions are those for which the order k is 
either an integer or a half-integer. Function values 
are quite difficult to compute for k a half-integer 
and x positive. 

Previously, interpolations from data compiled 
in tabular form were used to generate a compatible 
pair of Chebyshev approximations for k = 1/2. 


This allowed easy computation of F 1/2 ( x ) with a 
maximal relative error less than 5 x 10 -4 . 

Portions of the Loo Walsh arrays of rational 
Chebyshev approximations for k = -1/2, 1/2, 
and 3/2 are now available in tables. Maximal 
relative errors vary with the function and interval 
considered but generally range down to 10 -9 or 
less. 

Source: W. J. Cody and H. C. Thacher, Jr. 

Argonne National Laboratory 
(A RG- 10454) 

Circle 3 on Reader s Service Card . 


TABLE OF THE SOLUTIONS OF a TAN (ttx) = -b TAN Uttx) 


A reference table has been compiled which 
gives the first ten positive values of x satisfying the 
equation ofatan(7rx) - -b tan (a?rx) to four places 
for several values of a and b in the range 0.001 
< a < 1.0 and 0.001 < b < 1000. 

These solutions are used in the series solution 
of a problem involving transient diffusion in a 
composite slab. An eigenfunction technique is 
employed in this problem with the eigenvalues 
given by tan X^ = - M/a tan a X^, for k = 0,1,2,... . 

The transient response of two slabs of finite 


thickness in contact at a plane interface and 
initially containing nonequilibrium amounts of a 
transferable solute is of practical importance in 
the analysis of porous media in oil fields, and in 
the approach to equilibrium of a solute distribut- 
ing between two immiscible liquid layers. 

Source: P. Concus and D. R. Olander 
Lawrence Radiation Laboratory 
(LRL- 10007) 

Circle 4 on Reader's Service Card . 


REDUNDANT FORCE PROGRAM TO SOLVE VIBRATION RESPONSE PROBLEMS 


Missiles and space vehicles during static firing 
and flight experience dynamic vibrations which 
are usually sufficiently severe to degrade the struc- 
tural integrity. An analytical method that reliably 


predicts the response to a random environment 
has been developed and is available. 

This method is used in the analysis of redundant 
structures to obtain the dynamic response for 
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sinusoidal and random excitation. Random mo- 
tion of the base and multiple random loadings are 
considered. For example, the response to random 
motion of the base, an acoustic field, and both 
excitations simultaneously may be determined. 

Equations are developed for computing the 
response using generalized parameters. Cross 
power spectral densities are used to account for 
the correlation of the multiple random inputs. 
A lumped parameter idealization is employed in 
the method and structural influence coefficients, 
the natural frequencies and modes, and the dy- 
namic response of the structure are resolved. 

For the purpose of analysis, the required 
information consists of mass, damping, and 
stiffness matrices which may be obtained in a 
large variety of methods. If facilities are available, 
these parameters may be obtained experimentally. 

Once the transfer functions have been com- 


puted, the choice of the type of excitation to 
which the structure is subjected and the type of 
response to be computed is nearly unrestricted. 
The results include the accelerations, motions, 
and structural component loads for the system 
being analyzed. 

Versatility is the main advantage of this method. 
Its efficiency is demonstrated when used in con- 
junction with digital computers. The experienced 
user has the capability of increased accuracy by 
the use of more normal modes or increased ef- 
ficiency by using the most significant normal 
modes when computing the total response. 

Source: C. M. Fuller and D. N. Roudebush of 
McDonnell-Douglas Corp. 

under contract to 
Marshall Space Flight Center 
(MFS-12319) 

Circle 5 on Reader's Service Card. 


A WAVE TRANSMISSION METHOD OF ANALYSIS OF A TIMOSHENKO BEAM 


A better insight is gained into the problem of 
transverse bending of a beam having secondary 
dynamic effects of rotary inertia and shear com- 
pliance by the wave transmission approach. It 
provides a physically realistic representation of 
the structure and more accurate solutions than 
the standard-lumped model achieved through 
modal techniques. This is true since structural 
parameters are distributed and not lumped in 
actual conditions. 

The distributed approach results in an overall 
transmission matrix relating a set of system 
variables at one end of a beam to those at the other 
end, or may represent segments of the structure 


with variation in beam parameters. Also, two end- 
effect matrices and two propagation matrices are 
obtained through the solution of the Timoshenko 
beam equation in terms of the distributed system 
concepts of propagation and reflection, and 
through a transformation technique relating the 
local state variables to the characteristic variables 
of the beam. 

Source: W. C. Nowak of 
McDonnell Douglas Corp. 

under contract to 
Marshall Space Flight Center 
(MFS-12808) 

Circle 6 on Reader's Service Card. 


THE USE OF FOURIER TRANSFORMS IN THE SOLUTION OF BEAM 

AND PLATE PROBLEMS 


Fourier transforms can be used advantageously 
in the solution of static and dynamic problems of 
structural elements. The ordinary or the partial 
differential equation which describes the boundary 
value and/or initial value problem together with 
its prescribed boundary and initial conditions can 


be simplified by means of an integral transform 
and solved exactly, in closed form. This technique 
reduces a partial differential equation with n in- 
dependent variables to (n - 1) independent vari- 
ables and one transform parameter if the trans- 
form of the original differential equation is taken 
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with reference to one variable. By successive 
applications of the transform (n - 1) times with 
respect to the remaining (n - 1) transform vari- 
ables, one obtains an ordinary differential equa- 
tion with (n) transform parameters. Obtain the 
solution of the ordinary differential equation and 
use the transforms of the boundary conditions 
and then obtain the inverse. In the general trans- 
form approach (which includes Laplace, Mellin, 
Henkel and complex Fourier) the inverse is an 
integral and is not easy to solve. However, in finite 
Fourier Transforms, the inverse is a summation. 
However, sometimes it is more convenient to take 


n transforms with respect to n independent vari- 
ables. The various steps in the solution are shown 
in detail in the figure. 

Details of the general principles and some 
applications emphasizing Fourier finite sine trans- 
forms and Fourier finite cosine transforms only 
are available. 

Source: C. L. Amba-Rao of 
Wyle Laboratories 
under contract to 
Marshall Space Flight Center 
(MFS-12555) 

Circle 7 on Reader s Service Card . 


DYNAMIC BEHAVIOR OF CIRCULAR PLATES AND BEAMS 
WITH INTERNAL AND EXTERNAL DAMPING 


Nuclear reactors give rise to numerous prob- 
lems involving the dynamic behavior of elastic 
systems. Some of the more prominent concerns 
involve rapid startup or shutdown of reactors; 
control-rod scram; rapid opening and closing of 
pressure valves; and fluid-structure interaction 
in the core^and heat transfer system. In the case 
of mobile reactor and space vehicles, where weight 
is at a premium, complete system failure can re- 
sult if vibration problems receive inadequate 
attention. The development and use of some new 


versions of finite integral transforms in solving 
structural dynamic problems associated with 
nuclear reactor design are straightforward. 

A report is available which considers the 
dynamic response under various boundary condi- 
tions of circular plates including a rigidly clamped 
plate, simply supported plate, and clamped plate 
with an edge moment; the dynamic response of a 
cantilevered beam; and the dynamic interaction 
of a circular plate and cantilevered beam. 

The concept of mechanical impedance provides 
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a means for determining when it is possible to use 
an isolated elastic system in reactor design. For a 
single-degree-of-freedom system, the definition of 
impedance is the ratio of applied force to velocity 
for sinusoidal excitation. 

Use of the finite transform technique enables 
the solution of dynamic problems to be expressed 
in terms of normal modes. The solution appears as 
the sum of an infinite number of single-degree-of- 
freedom systems and facilitates physical inter- 
pretation and calculation. For lightly damped sys- 
tems, the first-mode term approximates the solu- 
tion to within 2% accuracy. 

The main advantage of the finite transform 
approach is that it handles static, forced vibration, 
aperiodic loads and time-dependent boundary 
conditions in a direct and simple manner, while 


avoiding the complex integration usually as- 
sociated with the Laplace transform technique. 

Use of the impedance concept can play an im- 
portant role in the analysis and synthesis of 
dynamic design problems and can simplify the 
solution of coupled elastic systems. 

The displacement (stress) is limited solely by 
the amount of damping in lightly damped systems. 
Inclusion of internal damping leads to the result 
that the stress is dependent upon strain rate as well 
as strain, even for systems remaining within the 
elastic limit. 

Source: G.Cinelli 
Argonne National Laboratory 
(ARG-10179) 

Circle 8 on Reader s Service Card. 


DYNAMIC VIBRATIONS OF THICK-WALLED ELASTIC ANISOTROPIC 
CYLINDERS AND SPHERES WITH INTERNAL DAMPING 


The dynamic response of thick elastic aniso- 
tropic cylinders and spheres is encountered in 
problems involving such diverse fields as pressure 
vessels, solid rocket propellants, geophysics, and 
ultrasonics. New integral transforms are used in 
determining the transient displacement and stress- 
es of material with transverse curvilinear ,isotropy 
and internal viscous damping for plane-strain 
motion of an infinitely long circular-cylindrical 
shell, torsional motion of a finite circular-cylin- 
drical shell, and radially symmetric motion of a 
spherical shell. Loads on the surfaces are taken as 
arbitrary functions of space (torsional case only) 
and time. 

The solution for the isotropic body is found by 
specifying certain values for the elastic constants. 
Cases of the solid body, thin shell, and cylindrical 


cavity in an infinite medium are obtained as limit- 
ing cases of the thick shell solution. An extended 
Weber transform is introduced which permits 
solution of infinite media problems in a more 
direct manner. 

By specializing the surface tractions to standard 
forms such as an impulse, step function, sinusoid, 
etc., the free, harmonic, and static motions are 
recovered. Resonance and mechanical impedance 
are derived and studied. Radially symmetric 
motion of a sphere is shown to be analogous to the 
plane-strain motion of the cylinder. 

Source: G. Cinelli 
Argonne National Laboratory 
(ARG-10178) 

Circle 9 on Reader’s Service Card. 


PREDICTION METHODS FOR FRICTION COEFFICIENTS FOR GASES 


Mathematical formulas are available for cor- 
relating both laminar and turbulent friction coef- 
ficients for gases, with large variations in the 
physical properties, flowing through smooth 
tubes. In recent years, the generally accepted rela- 
tions used to predict friction coefficients for 


laminar and turbulent flow have undergone 
extensive experimental investigation. Observed 
coefficients were found to be as much as three 
times the predicted value. 

Since 1967, empirical relations have been in 
use for laminar and turbulent flow wherein the 
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friction coefficient based on bulk temperature is 
a function of modified Reynolds number based 
on surface temperature. 

For modified surface Reynolds numbers less 
than 3000 (corresponding to laminar flow), local 
and average friction coefficients can be predicted 
within ±20 percent for 1< T s /T b < 4.1 by the 
expression: 

1 8 

2 = Re s 

where: 

f = friction coefficient 
T s = surface temperature 
T b = bulk temperature 
Re> = surface Reynolds number 
For modified surface Reynolds number of 3000 
and greater (corresponding to turbulent flow), 


local and average friction coefficients for both 
cooling and heating (0.35 < T s /T b < 7.4) can be 
predicted within ±10 percent by the expression: 

H °'° 007 + 05 

These relations have been used to correlate 
friction coefficients for hydrogen, helium, nitro- 
gen, carbon dioxide and air. 

These prediction methods apply to any of the 
problems of predicting friction pressure loss 
through a passage, from the most simple piece 
of tube to the complex cooling passages of a 
regeneratively-cooled rocket nozzle. 

Source: M. F. Taylor 
Lewis Research Center 
(LEW-10774) 

Circle 10 on Reader's Service Card . 


FREQUENCY DOMAIN ANALYSIS AND SYNTHESIS OF LUMPED PARAMETER 
SYSTEMS USING NONLINEAR LEAST-SQUARES TECHNIQUES 


A generalized technique for the design and 
analysis of linear systems in the frequency domain 
has been developed that greatly improves the 
design procedure for filter and compensation cir- 
cuits. The evaluation of lumped parametric system 
models has been found to be simplified and com- 
putationally advantageous in the frequency 
domain. 

The numerical Fourier transformation equa- 
tions involved in the frequency domain model are 
ideally suited for computer calculations. A 
nonlinear least-squares computer program has 
been developed which finds the least-squares best 
estimate for any number of parameters in an 
arbitrarily complicated model. 


Additional advantages of frequency domain 
modeling include: the facility with which ana- 
lytical frequency solutions can be obtained as 
compared with time solutions; there is no change 
in the form of the frequency solution as the values 
of the parameters are changed; and the limits of 
the integral of the squared error can be set to 
include only the region of interest. 

Source: J. R. Hays of 
The Boeing Co. 
under contract to 
Marshall Space Flight Center 
(MFS-15033) 

Circle 1 1 on Reader s Service Card . 


ESTIMATING RELIABILITY BY APPLICATION OF MATRIX REPRESENTATION 


In estimating the reliability of a manned space- 
craft for a predetermined mission, there are two 
probabilities to calculate: (1) mission success- 
completing the mission objectives with no loss of 
life, and (2) crew safety— returning the crew safely 
to earth with or without having accomplished the 


objectives of the mission. The large number of 
factors that must be considered when estimating 
these probabilities, such as the changing hardware 
configuration and the numerous components in 
a manned space vehicle, make computations 
complicated. 
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A new and simplif ying technique for performing 
this calculation has been developed, based upon 
matrix collapsing. This technique provides analytic 
expressions of mission success and crew safety 
for each subsystem, relating changes in subsystem 
reliability directly to mission success and crew 
safety. 

The reason for the development of the concepts 
of matrix representation and matrix collapsing is 
that ordinary Boolean techniques are cumbersome 
in cases where the configuration of a complicated 
subsystem is permitted to change with respect to 
time. These configuration changes may be due to 
changing functional or reliability requirements 
on the subsystem. 

The functional and reliability requirements of 
any given subsystem, “i”, for any phase, K,are 
represented in matrix form in order to facilitate 
deriving an analytic reliability function for the 
subsystem which is valid form time zero to the end 
of that particular phase. Even though there are 
three mutually exclusive states to consider — mission 
success,: abort enabling failure, and fatal failure— 
the only two matrices which need to be determined 
are the success matrix which contains all of the dis- 
tinct success paths through the subsystem, and the 
abort enabling failure matrix which contains all 
of the distinct abort enabling failure paths through 
the subsystem. The probability of a fatal failure 
at the end of any phase K is the probability of 
success at the end of phase K-l minus the prob- 
ability of crew safety to the end of phase K, with 


the total probability of a fatal failure being the 
sum of all the phase fatal-failure probabilities. 
The problem of determining the unconditional 
success matrix can be solved by utilizing the con- 
cept of matrix collapsing. 

Given two success matrices, A and B, where A 
is valid for the first mission phase and B is valid 
for the next mission phase in which the subsystem 
occurs, every row of A is intersected with every 
row of B and all redundant paths are eliminated. 
This process is continued for all subsequent mis- 
sion phases in which the subsystem occurs; this 
method is the same as intersecting the Boolean 
phase models, but is more simple. The abort 
enabling failure ground rules are applied to the 
unconditional success matrix to determine the 
unconditional abort enabling failure matrix. In 
order to include the aborts, the unconditional 
success matrix and the unconditional abort 
enabling failure matrix are collapsed with the 
subsystem functional and reliability requirements 
abort matrix. Once the unconditional matrices are 
determined, a reliability measure may be applied 
to them; the necessary functions are then derived 
in order to determine the probabilities of mission 
success and crew safety. 

Source: W. L. Austin of 
General Electric Co. 

under contract to 
NASA Headquarters 
(HQN- 10246) 

Circle 12 on Reader's Service Card. 


MATHEMATICAL DETERMINATION OF PARTICULATE CONTAMINATION 
LEVELS FOR SURFACE CLEANLINESS OF FLUID SYSTEMS 


A hypothesis that small particles exist as partic- 
ulate universes with a geometric mean size of 
approximately one micron applies in the paint, 
ceramic, and coal industries and serves as the 
basis for previous limitations on sizes and quanti- 
ties of particulate contamination permitted on 
cleaned surfaces of fluid systems. However, this 
hypothesis does not apply to complex systems 
which are unique and have their own requirements 
for definining what can be tolerated in fluids. A 
new method for the determination of particulate 
contamination levels for use in complex systems 


requires: (a) the definition of aparticle by a mathe- 
matical model; (b) a method for calculating the 
degree of contamination that can be tolerated; 
and (c) a method for estimating the probability 
that the contamination occurring on a surface 
will migrate with the fluid in the system. 

(a) Mathematical model: A projected diameter 
is defined as the diameter of a circle having the 
same area as the projected image of the particle 
when viewed through a microscope in the direc- 
tion perpendicular to the plane of greatest stabil- 
ity. The one-dimensional projected diameter of a 
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particle is related to a three-dimensional sphere by 
means of an equation employing a geometrically 
equivalent shape factor. The particle is sub- 
sequently calculated as an equivalent sphere by a 
conventional mensuration formula. 

(b) Contamination tolerance: Operational re- 
liability of a system is determined by the magni- 
tude of contamination deposited at critical clear- 
ances of mating parts and by the capability of 
these parts to continue to function. If a critical 
clearance is considered equivalent to a single pore, 
filter” media filtration mode equations can be used 
to derive an equation for correlating the volume 
of a critical clearance to the number of particles 
of a specific maximum size that can be tolerated 
to be deposited at that clearance. 

(c) Estimation of probability: Consider each 
particle an event in that class. The disposition of 
each particle (event) is mutually exclusive. Either 


the particle remains on the surface of the cleaned 
system or the fluid in the system removes it from 
that surface. Once entrained in the fluid system, 
the particle poses the threat of being deposited at 
a critical clearance. Whether the particle “stays” 
on the surface or “goes” with the fluid determines 
whether the event is favorable or unfavorable, 
respectively. Probability equations are used to 
calculate the probability of an unfavorable out- 
come for a single trial and the mean number of 
unfavorable events for “n” trials. Limitations can 
then be placed on the probability of an unfavor- 
able outcome. 

Source: Hayes International Corp. 

under contract to 
Kennedy Space Center 
(KSC- 10267) 

Circle 13 on Reader s Service Card . 


ACOUSTIC WAVE ANALYSIS 


A study has been made of the mechanism of 
generation of acoustic waves in a rotor/stator 
combination with a comparison of the relative 
strengths of waves initiated by different sources 
in a centrifugal pump. 

The velocity leaving a blade row is nonuniform 
due to both the gradients in the region designated 
as the potential flow and the viscous boundary 
layers along the blade surfaces, but is periodic 
in the mean, the period extending from blade to 
blade. This nonuniform velocity would appear to 
a downstream blade row in relative motion to the 
upstream blade row as a periodic unsteady veloc- 
ity field which produces unsteady pressures on 
the downstream blades. Acoustic waves are gen- 
erated by both the unsteady velocity field at the 
entrance of the blade row, and the unsteady blade 
pressure loads on the downstream blades due to 
the unsteady velocity field. 

A primary effort was made to indicate the 
relative importance of the acoustic waves gen- 
erated by these two related mechanisms. Before 
this could be done, the amplitude of the unsteady 
pressures on the downstream blade row had to be 
calculated. In the literature, these unsteady pres- 
sures have been investigated as a function of four 


effects which have been designated as the circu- 
lation, blade thickness, wake, and wake distortion 
effects. The analytical development of each is 
based primarily on the two-dimensional theory 
of the unsteady flow about a thin airfoil. Only 
the viscous wake effect was adopted for use in the 
current program, a computer program being writ- 
ten to calculate the unsteady pressures on a blade 
due to an approaching arbitrary wake velocity. 

The generation of waves into a medium sur- 
rounding a plate, by unsteady forces on the plate 
surface, is a coupled problem requiring solution 
of the equation of motion for transverse dis- 
placement of the plate, as well as solution of the 
wave equation in the medium. The two solutions 
are coupled by the boundary conditions which 
require continuity of both the normal velocity and 
pressure across the fluid-plate interface.” To 
estimate the order of magnitude of the wave 
amplitude generated by this coupled motion, a 
simpler example problem was considered which 
consisted of a plate of infinite extent and uniform 
thickness. The fluid medium above the plate was 
assumed to be coupled to the plate motion, but 
the fluid below the plate was not. The coupled 
equations were solved, yielding an expression for 
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the amplitude of the waves generated by an un- 
steady force on the plate. 

To estimate the order of magnitude of the 
amplitude of the waves generated by an unsteady 
velocity boundary condition at the inlet of a blade 
row (or duct), a second simple example was 
assumed which consisted of a straight, infinite, 
rectangular duct with a simple harmonic velocity 
at its inlet boundary, which velocity was uniform 
in the plane of the duct. The wave equation was 
solved to yield the generated wave amplitude. 

Using data for a standard centrifugal pump, the 
amplitudes of the oscillations generated by these 
two mechanisms were computed and compared. 
This comparison indicated that the amplitudes 


of the waves generated by the pressure loading on 
the blades (which in this case were volute tongues) 
were less by a factor of 10 4 to 10 5 than those gen- 
erated by the velocity boundary condition. Thus, 
the primary mechanism for generation of acoustic 
waves in the centrifugal pump, due to the rotor/ 
stator interaction, is that of an unsteady source 
at the entrance of the blade row as represented by 
the unsteady velocity field. 

Source: E. D. Jackson of 
North American Rockwell Corp. 

under contract to 
Marshall Space Flight Center 
(MFS-18076) 

Circle 14 on Reader s Service Card. 


ANALYTICAL DRAFTING CURVES PROVIDE EXACT EQUATIONS 

FOR PLOTTED DATA 


A rapid and accurate method has been devel- 
oped for obtaining explicit mathematical expres- 
sions for any numerical data that appear in the 
form of graphical plots, i.e., thermodynamic prop- 
erties of gases from a Mollier chart, or physical 
properties of materials. 

Such a method allows a considerable saving in 
time and effort when large quantities of tabulated 
data are needed in a particular computing pro- 
gram and there is insufficient storage space in the 
computer. An almost unlimited amount of data 
can be accurately represented by an exact mathe- 
matical equation and two transformation equa- 
tions. 

Given a set of data points that have been 
plotted in the coordinate system X t -Y t , one or 
more of the analytical drafting curves can be so 
oriented as to fit the plotted data to any desired 
degree of accuracy. In the drafting curves coor- 
dinate system X r -Y r , there results an exact mathe- 
matical expression or set of expressions for the 
plotted data. By use of the transformation equa- 
tions that involve translation and, in general, 
rotation, the plotted data can be represented in 
the primary coordinate system X r Y t . Therefore, 
the transformation equations (which are simple 
relations and are well known in mathematics) 
plus the mathematical equation of the drafting 


curve can be used to analytically represent all of 
the plotted data. 

It has been found from actual tests with second 
degree curves that the aforementioned method 
can easily provide agreement between the plotted 
data and the mathematical equations of better 
than 99% (where the plotted data are accurate). 
Experimental data, of course, would be within the 
accuracy with which a smooth curve represented 
the actual values. 

Tabulated data from solutions of high order 
equations often cannot be conveniently subjected 
to “curve fitting” by computer. Such data in 
graphical form, however, can be conveniently 
represented by the analytical drafting curve 
method. This method is most conveniently ap- 
plied when the plotted data are single valued in 
both coordinate systems. The use of analytical 
drafting curves has application to the teaching of 
analytic geometry, since two dimensional coor- 
dinate transformations and curves that result 
from equations with multiple roots can be graph- 
ically demonstrated in a unique manner. 

Source: R. B. Stewart 

Langley Research Center 
(LAR-90285) 

Circle 15 on Reader's Service Card. 
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NUMERICAL INTEGRATION OF ORDINARY DIFFERENTIAL EQUATIONS 

OF VARIOUS ORDERS 


Mathematical studies have resulted in tech- 
niques for the numerical integration of differential 
equations of various orders. Modified multistep 
predictor-corrector methods for general initial- 
value problems are discussed and new methods are 
introduced. Preliminary numerical evidence in- 
dicates that this process is faster and more ac- 
curate than similar methods applied to equivalent 
systems of first-order equations. 

The investigation was oriented towards dealing 
with higher-order equations directly, rather than 
converting them to a large system of first-order 
equations. It was found that the Nordsieck-type 
methods (N methods) were the most applicable to 
higher-order equations since the derivatives are 
directly available. In pursuing the relation between 
multistep methods and N methods, a large class 
of modified multistep methods became apparent. 
These k-step methods can have degree 2k, be 
strongly stable, and yet only use one derivative 
evaluation. 

An equivalence relation between the methods 
was determined. The Nordsieck formulation was 
found to be best for high-accuracy problems, but 
for low-accuracy problems either the N methods 
or versions of the new methods should be used. 
The conventional Adams multistep methods 
should almost never be used. 

The system of equations studied consisted of s 


equations in the s unknowns yj, their first pi deriv- 
atives, and the independent variable x. The most 
general case for which concrete results are given is 
when these equations can be solved for the high- 
est-order derivatives of each yi. 

For didactic reasons, the first sections of the 
study deal with the simplest case, s = p = 1, in 
order to prepare for discussions of the advantages 
of various formulations. Later, the N methods 
are applied directly to the above equation. A 
numerical comparison of the results of integrating 
the second-order equation describing J 16(X) 
directly and integrating the related pair of first- 
order equations is reported. Evidence indicates 
both a time and accuracy advantage in the direct 
approach. Some preliminary tests on the integra- 
tion of singular families are also reported. 

Finally, existence and convergence proofs are 
presented. The matrix formulation used provides 
a convergence proof equivalent to the usual proof, 
but considerably more compact and applicable 
to the system. 

Source: C. W. Gear 
University of Illinois 
under contract to 
Argonne National Laboratory 
(A RG- 10247) 

Circle 16 on Reader s Service Card . 


STRUCTURE OF THE ISOTROPIC TRANSPORT OPERATORS 
IN THREE INDEPENDENT SPACE VARIABLES 


The use of elementary functions for construc- 
tion of solutions for more than one space dimen- 
sion, if only for special' geometries, had been 
previously attempted. In this work, a decomposi- 
tion of the transport operator in three space 
dimensions is carried out. 

Based on the idea of separation of variables, 
a spectral theory for the three-dimensional, sta- 
tionary, isotropic transport operator is developed 
in a vector space of complex-valued Borel func- 
tions resulting in continuous sets of regular and 
generalized eigenfunctions. 

Because of the non-self-adjoint nature of this 


operator, the results could not be anticipated in- 
tuitively from the known decomposition of this 
operator in the special case of plane geometry. 
The results obtained indicate a promising nev) ap- 
proach to analytic solution of the linear transport 
equation in higher space dimensions. Examples 
are given for slab geometry with and without 
axial symmetry, spherical symmetry, and cylin- 
drical symmetries. 

This method is suggested by the principle of 
separation of variables, which is widely used for 
the decomposition of linear partial-differential 
equations of mathematical physics. The solution 
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space for the transport equation is mathematically 
defined as a set of complex valued Borel func- 
tions. For an attempt at completeness it was neces- 
sary to add to the set of regular eigenfunctions, 
and the corresponding regular expansion modes 
over a complex continuum, a set of generalized 
eigenfunctions and the corresponding generalized 
expansion modes over a complex continuum of 
even higher dimension. The regular eigenfunctions 
are characterized by a complex vector called a 
regular (or proper) eigenvector, A. These A 
eigenvectors, however, must satisfy the condition 
A* A = Aq, where assumes exactly one value for 
any given equation; therefore ±x 0 can be called 
eigenvalues. For yield of the generalized eigen- 
functions the improper integral must be evaluated. 

A basic result of this research is the general 
solution to the transport equation in a sense 
analogous to that of the theory of differential 
equations. Special representations are obtained 


by imposition of geometric restrictions on the 
general solution. These enable the numerical 
analyst to construct exact test problems for 
checking the performance of general multidimen- 
sional transport codes, and therefore can be most 
helpful in error-analysis. The calculation of 
singular integrals is not difficult if properly 
executed. The results of the work can be gener- 
alized to anisotropic scattering, multiregion, and 
multigroup problems. A similar theory can be 
conceived for the discrete ordinate approxima- 
tions; their spectra, modes, and general solutions 
can then be compared with the continuous case 
for a measure of fitness for the approximate solu- 
tions. 

Source; E. H. Bareiss and I. K. Abu-Shumays 
Argonne National Laboratory 
(ARG- 10448) 

Circle 17 on Reader's Service Card. 


ANALYSIS OF MAGNETICALLY-CONTROLLED PROCESSES 
IN PULSE-MODULATION SYSTEMS 


Simple analytic expressions have been estab- 
lished which are applicable to the design of a class 
of pulse modulators in which the modulating 
signal controls the reset level of flux in a nonlinear 
magnetic core. Analysis was carried out in the 
reset-coordinate plane, B r -T r , by defining two 
constraints on the reset process; one imposed by 
the magnetic core itself, and the other which 
includes the circuitry external to the core. Expres- 
sions derived for modulator response, modulation 
sensitivity, and limit-cycle stability were applied 
to both pulse-width and pulse-rate modulator 
designs. 

Basic character of modulator behavior was 
shown to depend on the relative magnitudes and 
algebraic signs of the circuit derivative dB r /dT r 
and the core partial derivative ^B r /^T r . Mod- 
ulator response was found to be oscillatory for 
designs in which (^Br/< rj, T r )/(dBr/dT r ) < 0, and is 
monotonic for (^B r /^T r )/dB r /dT r ) > 0. Mod- 
ulator stability depends only on the magnitude 
of (ciB r /^T r )/(dB r /dT r ). Discontinuities in the 
fundamental mode of operation or transitions to 


a subharmonic mode were shown to result from 
the onset of instability. The region of stable opera- 
tion is extended and the modulator response made 
very fast for a design using constant reset time. 

Results of the investigation show that mod- 
ulator performance parameters can vary over 
a wide range in the region of normal operation, 
with a significant deterioration in response occur- 
ring as the threshold of instability is approached. 
Pulse-width modulator designs with both finite 
and infinite circuit derivatives were analyzed for 
operation under constant-voltage reset conditions. 
Essentially optimum modulator performance was 
obtained for operation under fixed-reset-time 
conditions (infinite circuit derivative). The results 
of the generalized analysis were also applied to a 
pulse-rate modulator operating under constant- 
current reset conditions. Both an analytical de- 
scription of the magnetic core reset function and a 
graphic representation of an experimentally de- 
rived characteristic were used in the analyses, thus 
illustrating the flexibility of the generalized ap- 
proach. 
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Performance parameters of each of the above 
modulator circuits were measured and compared 
with the theoretical and calculated values. In each 
case good agreement was obtained and the results 
of the generalized analysis thereby substantiated. 


Source: Leo J. Veillette 
Goddard Space Flight Center 
(G SC- 10241) 

Circle 18 on Reader's Service Card. 


FINITE ELEMENT FORMULATION FOR LINEAR THERMOVISCOELASTIC MATERIALS 


The finite difference equations in time and finite 
element matrix equations in space have been 
developed for general linear thermoviscoelastic 
problems. 

There have been numerous applications of the 
finite element technique, mainly to elastic and 
plastic static problems and to some steady-state 
dynamic problems. However, the extension of this 
technique to viscoelastic problems without using 
the elastic-viscoelastic correspondence principle 
has been accomplished only in a few relatively 
simple cases. It has been concluded that the 
elastic-viscoelastic correspondence principle can 
be evoked for rather special cases, i.e., when the 
material properties are independent of thermal 
changes. Since the properties of most viscoelastic 
materials are highly temperature-sensitive, the 
development of a general program should be 
based on the solution of integral equations in real 
time rather than the correspondence principle. 

A brief statement of the thermoviscoelastic 


field equations is followed by the development of 
the finite different equations in time and then by 
the finite element formulation in space. The equa- 
tions are derived for a general three-dimensional 
body but are applicable with minor changes to 
one- and two-dimensional configurations. Some 
attention is given to the experimental determina- 
tion of material properties and their use in analyt- 
ical work. An expansion of the experimentally or 
analytically determined material property func- 
tions in terms of exponential series leads to 
recurrence matrix equations, eliminating the 
problem of calculating at each time-step the 
history of material response. 

Source: E. Heer and J. C. Chen of 
Caltech/JPL 
under contract to 
NASA Pasadena Office 
(NPO-1 1229) 

Circle 19 on Reader's Service Card . 


SOME NUMERICAL METHODS FOR INTEGRATE G SYSTEMS OF FIRST-ORDER 
ORDINARY DIFFERENTIAL EQUATIONS 


Numerical methods for integration of systems of 
coupled first-order ordinary differential equations 
have been studied and developed. The extrap- 
olation methods of Burlirsch-Stoer and Neville 
are explained, and the results of tests comparing 
these extrapolation methods with th//. Runge- 
Kutta and Adams-Moulton methods are given. 
There are some circumstances under which the 
extrapolation method may be preferred. Results 
of double-precision tests of the Bulirsch-Stoer 
method show the efficiency of various orders of 
this method as a function of desired accuracy. 

The speed and accuracy of the methods were 
compared by programming the methods on the 


same computer (CDC-3600) and applying each 
to the same sets of equations. The numbers of 
times that the functions to be integrated had to be 
evaluated, for a desired accuracy, were used as a 
computer-independent and program-independent 
measure of time. The accuracies were compared 
by use of each method to integrate from given 
initial conditions to some final point at which the 
relative error was computed. 

The results of cases run by the four methods 
tested indicate that the efficiency of the solution 
is dependent on proper matching of the order of 
the method of solution to the desired accuracy. 
For example, in tests on the negative-exponential 
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equation, for a relative error of 10 " ,6 , the Neville 
method with M = 4 takes twice as long as the 
method with M = 6 and about 2.5 times as long 
as the shortest time with M = 10. The Runge- 
Kutta fourth-order method and the Adams-Moul- 
ton fourth-order method take about 160 times as 
long as the Neville method with M = 10, or about 
250 times as long as the Bulirsch-Stoer method 
with M = 10. Even at relative errors of 10 " 9 , the 
Adams-Moulton and Runge-Kutta methods take 
about 9 times as long as the Neville and Bulirsch- 
Stoer methods with M = 6. At the other extreme, 
if a rough estimate is desired and errors up to 


10 " 2 can be tolerated, Neville’s method with 
M - 10 takes almost 7 times as long as the fourth- 
order Adams-Moulton methods, and the methods 
of lower order are preferable. 

The double-precision tests of the Bulirsch-Stoer 
method show it to be quite powerful for work of 
great accuracy; it seemed especially capable of 
coping with changing difficulty of solution of the 
equations presented to it. 

Source: N. W. Clark 
Argonne National Laboratory 
(ARG- 10308) 

Circle 20 on Reader's Service Card . 


CONTROLLABILITY OF DISTRIBUTED-PARAMETER SYSTEMS 


The controllability of distributed-parameter 
control systems (those which can be described by 
partial differential equations) has been studied, 
resulting in a general theory to include those con- 
trol systems which cannot be described by ordi- 
nary differential equations. 

The study summarizes techniques applicable to 
control systems problems. The eigenvalue-eigen- 
function expansion method for the solution of 
homogenous boundary value problems (BVP) 
is used. Problems in which the control appears 
at the boundary are treated by converting the 
nonhomogenous BVP to an equivalent homo- 
genous BVP by introducing generalized func- 
tions. The generalization of the concept of control- 
lability of finite dimensional systems to infinite 
dimensional systems is given. The pseudo-inverse 
of a linear operator is defined which is a generali- 
zation of that of a matrix for finite dimensional 


spaces. The pseudo-inverse is then used to obtain 
minimum energy control for distributed-param- 
eter systems. This generalization includes results 
for finite dimensional systems which are available. 
In the infinite dimensional problem, it is necessary 
to solve for the eigenvalues and eigenfunctions of 
an integral operator. The necessary and sufficient 
conditions for the states which are reachable when 
the control is required to satisfy a norm constraint 
are given. The conditions are obtained by an ap- 
plication of the moment problem to distributed- 
parameter systems. These results are then used 
to obtain conditions for complete controllability. 

Source: C. J. Herget of 
University of California 
under contract to 
Marshall Space Flight Center 
(MFS-19429) 

Circle 2 1 on Reader s Service Card. 


LINEAR SYSTEMS OF EQUATIONS SOLVED USING MATHEMATICAL ALGORITHMS 


A new mathematical algorithm has been devel- 
oped for the solution of linear systems of equa- 
tions, AX = B, which preserves the integer proper- 
ties of the coefficients. The algorithms presented 
can also be used for the efficient evaluation of 
determinates and their leading minors. 

The methods described are computationally 
more efficient than the corresponding single-step 
Gaussian elimination technique currently being 


used. In addition, the magnitude of the coefficient^ 
in the transformed matrices is minimized. The 
algorithms have universal applicability for solving 
linear equations and are particularly adaptable to 
digital computer techniques. 

Source: E. H. Bareiss 
Argonne National Laboratory 
(ARG-10146) 

Circle 22 on Reader s Service Card. 
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SHOCK AND VIBRATION RESPONSE OF MULTISTAGE STRUCTURE 


An analytical and experimental study has been 
made on the shock and vibration response of a 
multistage structure. Lumped-mass, continuous- 
beam, multimode, and matrix-iteration methods 
were used in the analytical phase of the study. 
Experimentally, a special technique was used in 
conjunction with a mechanical shaker. The study 
was made on the load paths, transmissibility, and 
attenuation properties along a longitudinal axis 
of a long, slender type of structure with increasing 
degree of complexity including ring frames, lon- 
gerons, bulkheads, propellant fluid, payload mass 
(spacecraft), as well as multistage structures. Free 
vibration characteristics were analytically solved 
by lump masses and continuous beam approaches. 
A modal approach was applied to determine shock 
responses to various phases of different forms and 
durations for a multideg ree-of-freedom system. 


In the experimental phase of the study, a special 
technique was employed to produce pulses of 
varying durations by a mechanical shaker. The test 
results verified the analytical predictions to a 
satisfactory degree. 

Longitudinal forced response to various pulses 
of different forms and durations can be calculated 
easily for slender or multistage structures, as in 
bridge support towers, and multistory building 
supports.* Solutions obtained are peak responses 
and complete shock wave propagation along the 
axial stations within the structure. 

Source: S. Y. Lee, S. S. Tang, and J. G. Liyeos of 
North American Rockwell Corp. 

under contract to 
Marshall Space Flight Center 
(MFS-14972) 

Circle 23 on Reader s Service Card . 


AN ORTHONORMALIZATION PROCEDURE FOR MULTIVARIABLE 
FUNCTION APPROXIMATION 


In many types of scientific and engineering 
problems, a table of two or more columns of data 
occurs and it is often desirable to present this 
data in a more useful form. The usual methods for 
performing this task are the many different tech- 
niques of multivariable function approximation 
such as the least-squares procedures that require 
appreciable time to compute the coefficients. 

Where a function of several variables is given 
numerically in tabular form, an orthonormaliza- 
tion technique allows an approximation of the 
numerical data to be determined in a convenient 
functional form. The method requires much less 


computational work than the usual least-squares 
technique, and allows more easily controlled ac- 
curacy. In this technique, the speed and accuracy 
of coefficient computation are much improved. 
Additionally, a very clear and useful physical 
interpretation of the procedure is available to aid 
in the choice of terms to be included in the ap- 
proximating formulas. 

Source: H. L. Ingram 
Marshall Space Flight Center 
(MFS-1313) 

Circle 24 on Reader s Service Card. 


MINIMUM PERMISSIBLE LEAKAGE RESISTANCE ESTABLISHED 
FOR INSTRUMENTATION SYSTEMS 


When an instrumentation system has been 
exposed to the elements, its leakage resistance to 
ground is appreciably affected to the detriment 
of system precision. Previously, such exposed 
systems have been dried out until leakage re- 
sistance to ground approaches infinity. This often 
far exceeds the precision/frequency requirements 


of a given system and thus results in a waste of 
time and money. 

Mathematical formulas are used to determine 
if, and to what extent, a given system should be 
dried out to restore minimum permissible leakage 
resistance. 

One example is a system that uses transducers 
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associated with a strain gage bridge network. The 
power supply minus lead of the bridge network 
is grounded at the test stand. Leakage from one 
of the bridge output leads to the stand ground 
unbalances the bridge. If the amount of leakage 
during a test changes sufficiently, the test data 
are invalid. Invalid data are readily identified by 
a significant zero shift and significant difference 
between post-test and pretest electrical unbalance. 

In a pressure measuring system that measures 
one megohm leakage to ground, system precision 
must not deteriorate by as much as 20%, i.e., 
0.0005 (0.2 x 0.0025). To determine whether the 
system should be dried out or the test run without 
this precaution, the following formula is em- 
ployed: 

( 106 ) 2 ft 

Ad _ 

” 106 + 0.8(35,000) 

0.0005 

AR = 17,55012 

where: AR = Change in leakage resistance 
during test. 

This example shows that if the leakage resistance 
increased or decreased 17,550 ohms from the 
initial value of 1 megohm during test, system 
precision would deteriorate 0.05%. Because 17,550 
ohms represents only a 1.755% change, it is likely 
that the leakage resistance would change by this 
amount. In this case, the system would be dried 
out before running the test. 

Assuming the system dried out so that leakage 
to ground measures 10 megohms, the question of 


whether to dry out the system further or run the 
test is resolved by the following: 

ar im 2 n_ 

1Q7 + 0.8(35,000) 

0.0005 

AR = 1.5M12 or 15% change. 

In this case the test would be run without further 
drying since it is unlikely the leakage resistance 
would change as much as 15% during test. 

When it is assumed the leakage resistance will 
not change more than 20% during a test, and that 
a deterioration in precision of 20% of the pre- 
cision classification of the system can be tolerated, 
the following formula determines the minimum 
permissible leakage resistance: 

R 

L N P 

where: R L = Minimum permissible leakage re- 
sistance; R c = Resistance of electrical unbalance 
resistor for 80% deflection; N = Number of elec- 
trical unbalance resistors; P = System precision 
classification (0.0025, 0.005, 0.01, 0.02, etc.). 

Formulas may be derived to be used for an 
indeterminate number of instrumentation systems 
that are exposed to moisture penetration. 

Source: J. L. Perrin of 
North American Rockwell Corp. 

under contract to 
Marshall Space Flight Center 
(MFS-848) 

Circle 25 on Reader s Service Card. 


Section 2. Computer Programs for Industrial Applications 


COMPUTER PROGRAM (PI-GAS) CALCULATES P-0 AND P-1 TRANSFER 
MATRICES FOR NEUTRON MODERATION IN MONATOMIC GAS 


To generate multigroup thermal neutron scat- 
tering cross sections for use in neutron transport 
theory programs, it is necessary to have the group- 
to-group elastic transfer kernels. Below 1 eV, ex- 
perimental crystalline transfer kernels may be 
used, at about 1 eV, it is necessary to analytically 
generate the elastic transfer kernels. The problem 
is to generate the P-0 and P- 1 transfer matrices for 
neutron moderation in a monatomic gas. 


The equations programmed in the PI -GAS 
FORTRAN IV Program are essentially those 
derived by Clendenin and others. These equations 
for the elastic scattering kernels are based on the 
following conditions: (1) there is isotropic scat- 
tering in the center-of-mass coordinate system; 
(2) the scattering cross section is constant, i.e., 
independent of the velocities; and (3) the target 
nuclear velocities satisfy a Maxwellian distribu- 
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tion. The monatomic gas model has the advantage 
that the differential cross section (matrix elements) 
may be expressed analytically. 

To calculate the elastic scattering transfer 
matrices for a neutron in a solid moderator, one of 
three models may be used, depending on the 
energy of the neutron and the temperature of the 
moderator. When the neutron energy is greater 
than a few eV, the scattering process depends only 
weakly on the facts that the target nuclei are 
bound in the lattice and the target nuclei have 
thermal motion. Thus, the moderator nuclei may 
be treated as free and stationary. In this model 
the neutron energy cannot increase as a result of 
a collision. For lower neutron energies (the 
thermal range) “upscatter” can occur due to the 
thermal motion of the moderator nuclei. For neu- 
tron energies less than about 1 eV, the interaction 
of the neutron range (energies above about 1 eV), 
the gas model described here is applicable. In the 
gas model, no account is taken of the fact that the 
scattering nuclei are bound in ordered fashion in 
the lattice; however, the thermal motion of the 
target nuclei is considered. 

When neutrons are moderating in a monatomic 
gas such as helium, the results of PI-GAS are 
directly applicable. In this case the thermody- 
namic temperature of the gas characterizes the 


Maxwellian distribution. However, for a molecu- 
lar gas such as hydrogen, the results in general do 
not apply since the internal degrees of freedom 
(vibrations and rotations) are not treated in the 
gas model. 

Quantities are printed out in the program that 
provide a check of the calculation to ascertain 
whether or not specific conditions are satisfied by 
the transfer matrices. 

This program is written in FORTRAN IV and 
MAP for the IBM 7090/7094 computer. 

One of the input quantities to the PI-GAS pro- 
gram is the temperature of the gas. Lamb first 
showed (for the resonance absorption of neutrons) 
that a Maxwellian distribution may be used for the 
nuclear velocities in a Debye solid. Nelkin and 
Parks showed this was also true for slow neutron 
scattering. However, the Maxwellian distribution 
is not characterized by the thermodynamic tem- 
perature T 1 of the moderator, but an effective 
temperature T 1 directly. Actually, there are addi- 
tional methods for obtaining T 1 . 

Source: G. Gibson and G. Collier of 
Westinghouse Astronuclear Laboratory 
under contract to 
Space Nuclear Propulsion Office 
(NUC-10141) 

Circle 26 on Reader s Service Card . 


MOP (MATRIX OPERATION PROGRAMS SYSTEM) 


The programming system consists of a set of 
FORTRAN IV subroutines which are related 
through a small common allocation. The system 
accomplishes all matrix algebra operations plus 
related input-output and housekeeping details. 
It was’ codedfor the IBM 7094 DCOS, but the 
document contains specific instructions for con- 
version to other computers. It is a convenient and 
complete programming base for jobs using 
substantial matrix algebra manipulations, and in 
no way prevents the user from calling other sub- 
routines. 

The system can be used by writing a main pro- 
gram which calls the MOP system subroutines. 
After the system is initialized, all matrix algebra 
operations can be performed. The casual program- 
mer can operate in an automatic nominal mode; 


the professional programmer has access to all sys- 
tem control flags and may exercise subtle options 
within the system. The naming of calling se- 
quences and operations is done according to an 
easily mastered mnemonic scheme. 

A simple overlay scheme allows loading of the 
entire system in 6000 words of core, which, with 
the system subroutines and a substantial main 
program, leaves sufficient space for two 50 x 50 
square matrices with four symmetric operators 
overlaid for flexibility. The sophisticated house- 
keeping and internal conversions of the system 
may be loaded if desired, but high efficiency is not 
claimed for programs desiring only a few calls. 

The above applies to double precision word 
length, but the system can be rapidly converted 
to single-precision if desired. 
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This system is written in FORTRAN IV for use 
on the IBM 7094 computer. Timing and program 
accuracy performance is good, and complete de- 
bugging programs are included which automati- 
cally test the system. The checkout has been 
exhaustive, and the document is quite complete. 


Source: P. M. Muller 
NASA Pasadena Office 
(NPO- 10429) 

Circle 27 on Reader's Service Card. 


COMPUTER PROGRAM FOR MASS OPTIONAL SOLUTIONS OF SOME 
END-POINT TRAJECTORY PROBLEMS 


The calculation of trajectories is required for a 
variety of orbital problems such as orbit transfer, 
rendezvous, lunar transfer, and lunar launch. In 
all cases optimization of trajectories for minimal 
propellant consumption is a prime concern. This 
program incorporates the technique of calculating 
coast arcs into a three-dimensional fixed end-point 
steepest ascent computer program. There is no 
restriction on the magnitude of thrust or on initial 
or final orbit characteristics. The initial condi- 
tions and desired terminal conditions of a transfer 
trajectory can be specified in conventional orbital 
elements or in spherical coordinates. 

The equations of motion are evaluated in spher- 
ical coordinates whatever the form of the input. 
Hence, the program actually calculates a tra- 
jectory between any two points in space defined 
by initial and final position vectors. The cor- 
responding end-point velocities are equally 
arbitrary. Problems can be handled which, on the 
surface, do not appear to fall within the orbit 
transfer category. 

The problem of thrusting flight in a vacuum in 
the presence of two-body forces only is formulated 


in the calculus of variations and solved by the 
method of steepest ascent. Optimal thrust direc- 
tions and thrust durations are found for a variety 
of orbital transfer problems in two or three di- 
mensions. Constant thrust and specific impulse 
are assumed. 

The program can be used on either the IBM 
7094 or SRU 1107 and has been written in 
FORTRAN IV language. 

Aerospace, missiles, ballistics, physics, and 
geology have some general similarity in flight 
problems, either of projectiles or particles, to 
which this program is applicable. Also this tech- 
nique could be of value in other isolation tech- 
niques such as the calculus of variation or New- 
ton-Raphson. 

Source: M. S. O’Mahony, A. G. Bennett, 
and C. D. Eshridge of 
The Boeing Company 
under contract to 
Marshall Space Flight Center 
(MFS-12976) 

Circle 28 on Reader's Service Card. 


HICOV (NEWTON-RAPHSON CALCULUS OF VARIATION 
WITH AUTOMATIC TRANSVERSALITIES) 


In generating trajectories that are optimum with 
respect to payload placed in an earth orbit, there 
is a need for analytically forming the transversality 
equations and their respective partial derivatives. 

A computer program generates the desired 
trajectories through the use of a subroutine 
package which produces the terminal and trans- 
versality conditions and their partial derivatives. 

The Wierstrass condition defines the control 


variables, or steering angles, which appear in the 
differential equations of motion, while the 
transversality conditions are used along with the 
physical cutoff conditions to determine the nine 
boundary conditions. The partial derivatives are 
used to solve for the desired boundary conditions 
and then determine the gain matrices for the steer- 
ing angles. The gain matrices are used to write 
the guidance functions which reduce to linear 
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polynomials in the state variables. The state vari- 
ables can be stored onboard the vehicle for guid- 
ance purposes. 

The transversality generator program along with 
theFORMAC Subroutine package, TAC, generates 
the terminal and transversality conditions and 
their partial derivatives. The package uses a 
formula manipulation language which is an ex- 
tension of FORTRAN. It allows the programmer 
to arrive at a single analytic solution, following 
it by multiple numeric evaluations. 

The terminal conditions are input into the pro- 
gram and a FORTRAN IV source deck for the 
calculation of the solutions is generated on tape. 
Included in the tape source deck are a deck name, 
calling sequence, type statements, dimension 


statements, and analytic expressions as desired. 
The tape can then be punched on cards and used 
directly with a program or called in by a program 
through the use of $ IEDIT. 

Using this subroutine package, the programmer 
is able to automatically arrange his analytic solu- 
tions) in form suitable for numeric evaluation. 

This program is written in FORTRAN IV 
and FORMAC for the IBM 7094 computer. 

Source: T. J. Heintschel of 
General Electric Company 
under contract to 
Marshall Space Flight Center 
(MFS- 14468) 

Circle 29 on Reader s Service Card. 


LEAST-SQUARES SMOOTHING (LSSMG) PROGRAM 


This program was developed so that the smallest 
computation time could be used to obtain a 
“smoothed” value for 21, 41, 61, or 81 points of 
raw data. For a slightly longer time, any odd 
numbers of points < 81 may be used. 

For an odd-numbered set of ordered pairs (x, y), 
which should represent some approximation to a 
function y =f (x) which is at most second degree, 
LSSMG will automatically provide the medium 
pair (x m , y m ) and will also provide the derivative 
dy/dx at the point x =x m , computing so much of 
the function y = ai + a 2 x -j- a 3 x 2 

as is required to give ai and a 2 
( -dy/dx when x m - 0) 


Since the primary purpose is speed, the function 
is not computed, nor is the array of smoothed 
points computed, but these may be obtained by 
the calling program by means of additional FOR- 
TRAN coding. An error routine is included for 
n > 81 or n even. 

Language: GMAP Assembly Language 

Machine Requirements: GE-625 

Source: Wallops Station 
(WLP- 10031) 

Circle 30 on Reader s Service Card. 


CALCULATION OF THE VOLUME OF CAVITY ON AN INDUCER BLADE, 
ASSUMING A FINITE LEADING-EDGE CAVITY 


This program calculates cavity boundaries, such 
as the shape and cavity volume, for a cavitating 
inducer blade by evaluating the Free Streamline- 
Wake Theory. The complex variable potential- 
flow solution is modified so that the branch cuts 


needed for evaluation of the solution exactly match 
the branch cuts defining the Fortran H complex 
logarithm routine CLOG. The ratio of discharge 
velocity to the inlet velocity is expressed as a 
function of cavitation number K and blade geom- 
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etry, so that the solution may be evaluated for a 
series of K values and flow and blade angles. The 
volume of the cavity is found by integration of 
the net area between cavity boundary and blade- 
suction surface. 

Output consists of plots of the free streamlines. 
Also, the value of cavity area, metal area, net 
cavity length/ spacing ratio, and cavity height/ 
spacing are tabulated. 


Language: FORTRAN H 

Machine Requirements: IBM 360, Release II 

Source: North American Rockwell Corp. 

under contract to 
Marshall Space Flight Center 
(MFS-18093) 

Circle 3 1 on Reader's Service Card. 


BELLOWS CALCULATION PROGRAM 


This program employs empirical and analytical 
derived design equations on various metal bellows 
of different sizes in order to calculate various 
properties of bellows used in ducting systems. 
Arithmetic operations are performed in double 
precision. Calculations are restricted to four single 
bellows movements and two double bellows move- 
ments. One subroutine and one data deck are 
required with the main program. 

The main program and the subroutine calculate 
bellows spring rates, bulging, bending, and hoop 
stresses. Cycle life is calculated by data deck. With 
known bellows dimensions and type of movements 


supplied as data, main program and subroutine 
calculate spring rate, actuating force, squirming 
pressure, stress, bellows weight, resonant fre- 
quency, fatigue life, and convolution clearing. 

Language: FORTRAN H 

Machine Requirements: IBM 360, Release 11 
Source: North American Rockwell Corp. 

under contract to 
Marshall Space Flight Center 
(MFS-12641) 

Circle 32 on Reader s Service Card. 


CALCULATION OF ISOTHERMAL, TURBULENT-JET MIXING OF TWO GASES 


This program solves a simplified model of the 
turbulent mixing that occurs between a central 
fuel jet and a surrounding, faster-moving coaxial 
stream of propellant. The von Mises transforma- 
tion was used to convert the axisymmetric forms 
of the isothermal boundary layer momentum and 
diffusion equations to forms amenable to numeri- 
cal solution. The effects of confining walls were 
not considered. The program can solve problems 
in which the initial velocities and densities of the 
two streams differ greatly, by using expressions 
for eddy viscosity that vary radially as well as 
axially. 

The input to the program consists of the initial 
ratio of coaxial -stream velocity to jet velocity, the 
mass fraction of component one in the initial jet 
and in the initial coaxial stream, the ratio of 
molecular weight of component two to the molec- 
ular weight of component one, the constants in 


the eddy viscosity formulation, the turbulent 
Schmidt number, the ratio of reactor diameter to 
jet radius, the constants in the reference density 
formulation, and the axial positions at which out- 
put is desired. 

At each axial position, the program prints out 
the axial position, the eddy viscosity, the product 
of density and eddy viscosity, and (for axial veloc- 
ity, mass fraction, and mole fraction) the quantity: 

(Centerline value-Coaxial-stream value) 

(1 -Coaxial-stream value). 

The radial variations with stream function, also 
converted to radial position, and the ratio of radial 
position to half radius for axial velocity, mass 
fraction, and mole fraction are provided in the 
following form: 

(Local value-Ambient stream value) 
(Centerline value-Ambient stream value) 
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Also listed are momentum flux normalized with 
the square of the centerline velocity, mass flux 
normalized with the product of centerline velocity 
and mass fraction, and eddy viscosity and the 
product of density and eddy viscosity divided by 
their centerline values. Both sides of the centerline 
compatibility condition are printed next, followed 
by the values of velocity and density at the largest 
stream function position in the calculation. Finally, 


the “line-of-sight” concentration, the dimension- 
less mass of component one, and the ratio of mass 
of component one-to-initial mass are listed. 

Language: FORTRAN IV 
Machine Requirements: IBM 7094 

Source: Lewis Research Center 
(LEW- 10579) 

Circle 33 on Reader's Service Card. 
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